
	/********************************************************************/
	/*                                                                  */
	/*                                                                  */
	/*              Copyright (C), 1993  SEIKO EPSON Corp.              */
	/*              ALL RIGHTS RESERVED                                 */
	/*                                                                  */
	/*                                                                  */
	/*   file name : modf.c                                             */
	/*                                                                  */
	/*   Revision history                                               */
	/*       93/08/31    first release                 K.Shinozuka      */
	/*       93/12/17    change include file           K.Shinozuka      */
	/*       97/03/17     modified	    		   T.Katahira	    */
	/*			change function declaration		    */
	/*					to the ANSI style	    */
	/*                                                                  */
	/********************************************************************/

#ifdef UNIX
#include "math.h"
#include "smcvals.h"

#else
#include <math.h>
#include <smcvals.h>
#endif

double
modf( double x, double * nptr )

	/********************************************************************/
	/*                                                                  */
	/*  The modf function devides a floating point number into a        */
	/*  fraction and an integer parts, each of which has the same sign  */
	/*  as the argument.                                                */
	/*  The integer part is stored as a double in the object pointed    */
	/*  to by nptr.                                                     */
	/*                                                                  */
	/*  RETURN    normal: the signed fractional part of x               */
	/*            error : NaN (input error)                             */
	/*                                                                  */
	/********************************************************************/

/* ALGORITHM
 * 1.get the fraction part as xfrac and exponent part as xexp from the argument
 * 2.clear fractional bits, if x can be separate into the fraction and the integer
 * 3.get the fractional number by subtracting the integer from x 
 */

/*  the architecture of double floating point
 *
 *   0 1          11                  31 32                               63 bit
 *   -----------------------------------------------------------------------
 *  | |   exponent  |                    fraction                           |
 *   -----------------------------------------------------------------------
 *
 *  |               |                   |                                   |
 *  |    12 bits          20 bits       |             32 bits               |
 *  |            lower word             |            higher word            |
 *
 *         bit    0         sign bit         (  1 bit  )
 *              1 - 11      exponent part    ( 11 bits )
 *             12 - 63      fraction part    ( 52 bits )
 *
 */
 
	/* double         x;           input value                              */
	/* double         *nptr;       the pointer, the integer to be stored    */
	{
	 DBL_LNG       ckarg;       /* exchange the type to check the argument  */
	 DBL_LNG       ipart;       /* integer stored in nptr                   */
	 DBL_LNG       xfrac;       /* fraction part of x                       */
	 DBL_LNG       ifrac;       /* fraction part of integer                 */
	 DBL_LNG       ffrac;       /* fraction part of fraction                */
	 register long xexp;        /* exponential part of x                    */
	 register long iexp;        /* exponential part of integer              */
	 register long fexp;        /* exponential part of fraction             */
	 register int  sign;        /* sign of x, 1 ( x>0 ) or -1 ( x<0 )       */

/* check the argument */	
	 *nptr = NAN._D;                          /* domein error */
	 ckarg._D = x;
	 CHECK_ARG( &ckarg );                     /* check NaN */
	                                          /* error: error code 33 domain error */
 
	 if ( ( ckarg.st._LL | NEG_ZERO_L ) == N_INF.st._LL ) {   /* Inf or -Inf */
	    errno = ERANGE;                       /* 34: range error */
	    *nptr = x;
	    return( 0.0 );
	 }

	 if ( x < 0.0 ) {                         /* check the sign of x */
	    x = - x;
	    sign = -1;
	 } else if ( x == 0.0 ) {
	    *nptr = 0.0;
	    return ( 0.0 );
	 } else {
	    sign = 1;
	 }

/* separate x into fraction & exponent part */
	 xfrac._D = x;                
	 xexp = (( xfrac.st._LL & MASK_EXP ) >> L_FRAC );    /* get exponent part  */
	                                          /* MASK_EXP = 0x7ff00000 */
	 xexp = ( xexp - ((unsigned long)NUM1_L >> L_FRAC) );    /* get off offset */
	                                          /* NUM1_L = 0x3ff00000 */
	 xfrac.st._LL = ( xfrac.st._LL &  MASK_FRAC );    /* get fraction part  */
	                                          /* MASK_FRAC = 0x000fffff */
	 if ( xexp >= N_FRAC ) {                  /* no fractional bits  */
	    ipart._D = x;                         /* set integer */
	    x = 0.0;                              /* set fraction */
	 } else if ( xexp < 0 ) {                 /* x < 1 : all fractional bits */
	    ipart._D = 0.0;                       /* set integer */
	 } else {
	    ipart._D = x;
	    if ( xexp <= L_FRAC ) {               /* integral bits are only in lower word */
	       ipart.st._LL = ( ipart.st._LL & ~( (unsigned long)MASK_FRAC >> xexp ));
	                                          /* clrear fractional bits */ 
	       ipart.st._LH = 0;
	    } else {
	       ipart.st._LH = ( ipart.st._LH & ~( ~(unsigned long)0 >> ( xexp - L_FRAC )));
	    }
	    x = x - ipart._D;                     /* get fraction */
	 }
	 if ( ipart._D == 0.0 ) {                 /* -0.0 is not allowed */
	    *nptr = 0.0;
	 } else {
	    *nptr = (double)sign * ipart._D;
	 }
	 if ( x == 0.0 ) {                        /* -0.0 is not allowed */
	    return( x );
	 }
	 return( (double)sign * x );
	}
